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Abstract 

Using a three-dimensional mean-field model we study one-dimensional dipolar Bose-Einstein condensate (BEC) soli- 
tons on a weak two-dimensional (2D) square and triangular optical lattice (OL) potentials placed perpendicular to the 
polarization direction. The stabilization against collapse and expansion of these solitons for a fixed dipolar interaction 
and a fixed number of atoms is possible for short-range atomic interaction lying between two critical limits. The soli- 
tons collapse below the lower limit and escapes to infinity above the upper limit. One can also stabilize identical tiny 
BEC solitons arranged on the 2D square OL sites forming a stable 2D array of interacting droplets when the OL sites 
are filled with a filling factor of 1/2 or less. Such an array is unstable when the filling factor is made more than 1/2 by 
occupying two adjacent sites of OL. These stable 2D arrays of dipolar superfluid BEC solitons are quite similar to the 
recently studied dipolar Mott insulator states on 2D lattice in the Bose-Hubbard model by Capogrosso-Sansone et al 
[B. Capogrosso-Sansone, C. Trefzger, M. Lewenstein, P. Zoller, G. Pupillo, Phys. Rev. Lett. 104 (2010) 125301]. 
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1. Introduction 

After the experimental observation of Bose-Einstein condensate (BEC) Q, there has been great interest in the 
problem of stabilization of BEC on a periodic optical lattice (OL) potential. The OL potential is generated in a 
laboratory by a standing wave polarized laser beam [2]. The resulting periodic potential simulates the potential seen 
by an electron in a solid |3|. As one of the interests in studying BEC droplets of small number of atoms on OL 
is to generate a stable array of BEC droplets by occupying each OL site with one tiny droplet so that an array of 
matter wave is formed as in condensed-matter physics (4). Unlike in condensed-matter physics, these BEC droplets 
are completely pure. Such a pure array of matter- wave simulating a quantum solid can be created and studied in a 
laboratory to model many quantum-mechanical condensed-matter phenomena. 

The above study to model condensed-matter phenomena with matter wave has obtained new impetus after the 
observation of dipolar BEC of 52 Cr [6) and 164 Dy (7) with a large dipolar interaction. The dipolar interaction is 
anisotropic and of long range in contrast to the short-range isotropic atomic interaction. Because of the anisotropic 
long-range interaction, the conditions of stability of a dipolar BEC soliton follow a distinct trend from that of a 
nondipolar BEC soliton [ 8 ]. Also, dipolar atoms have permanent magnetic moment | 5 ] and if polarized by an external 
magnetic field, an array of tiny BEC droplets with magnetic moment can simulate the problem of generation of 
magnetism in solids from individual atomic or molecular magnetic dipoles. 

Here, using a three-dimensional (3D) mean-field model, we study one-dimensional (ID) dipolar BEC solitons, 
free to move along the polarization direction z, on a two-dimensional (2D) square and triangular periodic OL in the 
x-y plane. In both cases, for a fixed dipolar interaction and fixed number of atoms, the solitons are stable for atomic 
short-range interaction (scattering length) between two critical limiting values. Below the lower limit the dipolar 
and short-range interactions lead to too much attraction and the soliton collapses and above the upper limit the net 
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attraction is too weak and the soliton escapes to infinity. The spreading of the BEC soliton in the 2D OL in the x-y 
plane is stopped by the attractive dipolar interaction along polarization direction z - the system lowers its energy by 
being long in the z direction and thin in the x-y plane. 

First, we consider a dipolar BEC soliton on a 2D square or triangular OL using the numerical and Lagrangian 
variational analysis of the mean-field Gross-Pitaevskii (GP) equation. The variational results are found to be in good 
agreement with numerical results. To demonstrate the stability, we perform a linear stability analysis |9] using the 
variational solution and calculate the normal-mode frequencies ifTOl . 

In an attempt to generate an array of dipolar BEC solitons on a 2D square OL, each with a small number of 
atoms, we find that interesting stable periodic structure can be formed for filling factors of 1/2, 1/3, and 1/4 termed 
checkerboard, stripe, and star configurations, respectively. Similar stable structures of ultra-cold dipolar atoms were 
first obtained as Mott insulator states in a study of dipolar atoms on 2D OL by solving the corresponding field- 
theoretic 2D Hubbard model numerically by Monte Carlo technique fTfl . This suggests that such stable structures are 
a consequence of the typical repulsive dipolar interaction in the x-y plane. 

There have been studies of 2D dipolar BEC solitons, free to move in the x-z or x-y plane with harmonic traps 
along y fT2l or z lfl3ll axis, respectively, and of ID dipolar BEC soliton under transverse harmonic trap [14]. The 
present ID dipolar BEC soliton confined by only a weak 2D OL in the x-y plane is distinct. The BEC soliton of 
the previous studies [Q21 [T3J [14) will essentially have a Gaussian density distribution along the infinite trap direction, 
whereas the present BEC soliton will have an exponential density distribution due to weak finite OL traps in these 
directions. Similar nondipolar BEC solitons in a lower dimensional OL have also been studied fT5l . 



2. Analytical formulation 

We consider a ID dipolar BEC soliton of N atoms, each of mass m, using the GP equation: 

V 2 



<90(r, t) 



- y + v ol + AO + uddifldd, N) 



(r,0, (1) 



with the bulk chemical potential jd(a,N) = 4nan, n = A/]0| 2 , with a the atomic scattering length, n the density, 
and r = {x, y,z] = {p,z}. The dipolar bulk chemical potential is fidd(add,N) = N f Uddif - r)|0(r r , t)\ 2 dr f , C/^(R) = 
3a^(l - 3 cos 2 6)/R 3 the dipolar interaction potential, R = r - r', normalization J (p(r) 2 dr = 1, 6 the angle between 
R and polarization direction z, Vql = _y ocos(2x) - Vocos(2v) for the square OL, and =-VoCOs(2x) - Vocos(x + 
V3v) - Vo cos(x - V3y) for the triangular OL, with Vb the strength of the OL. This 2D triangular OL consists of three 
OL in the x-y plane at mutual angles ofn/3. The length add = j^o/2 2 m/(l27ih 2 ) is the strength of dipolar interaction, jl 
the magnetic dipole moment of an atom, and yUo the permeability of free space. In Eq. ([I]), length is measured in units 
of /o = A/2n, taken here as 1 /mi, and time t in units of to = mftjh, where A is the OL wave length. Energy E and Vo 
are expressed in units of 2E R = h 2 /(mA 2 ), where E R is the recoil energy of one atom of mass m absorbing one lattice 
photon of wave length A. 

First we consider a ID BEC soliton on the OL V 2 ^. The Lagrangian density of Eq. |TJ) is fT6l 

L = l - (#* - + i|V0| 2 + 2naN\cf>\ 4 + V 2 £|0| 2 + ^N\cf>\ 2 J U dd (r - r')|0(r')| 2 </r'. (2) 

For a variational study we use the Gaussian ansatz [16]: 0(r, t) = exp(-p 2 /2w 2 - z 2 /2w 2 +/yp 2 + i/3z 2 )/(w p Vw z 7r 3/4 ) 
where w p and w z are time-dependent widths and y and ft are time-dependent chirps. Because we consider identical 
strengths of OL along x and y directions, an axially-symmetric Gaussian profile for the density is a good approximation 
to the actual density. The effective Lagrangian L (per particle) is 
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with kinetic, trap, and interaction energies given, respectively, by E^ n = -[1 /(2w 2 )+l/(4w 2 )], E tmv = -2Vo exp(-w 2 ), 
for the square OL and = -3Vq exp(-w 2 ), for the triangular OL, E mi = N[a - addj '(*)]/( V27rw 2 w z ), where /(k) = 
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Figure 1: (Color online) (a) Numerical 3D contour of a nondipolar BEC soliton of 1000 atoms with the scattering length tuned to a = -4.5<zo 
(g = AnaN = -3,gdd = 3addN = 0). The same for a dipolar BEC soliton of 1260 52 Cr atoms with scattering length a tuned to (b) a = 
(g = 0, gdd = 3) and (c) a = \2clq (g = l,gdd = 3). (d) The same for a dipolar BEC soliton of 4200 52 Cr atoms with scattering length a tuned to 
a = 10.75<2o (g = 30, gdd = 10). The square OL potential is Vq£ = -2 cos(2x) - 2 cos(2_y). The density |0(r)| 2 on the contour is 0.001. 



[1 + 2k 2 - 3k 2 d(K)]/(l - k 2 \ d{K) = atanh Vl - k 2 / VI - k 2 , k - w p /w z . The Euler-Lagrange equations for parameters 
w p ,w z ,y,/3 can be used to obtain the following equations of the widths for the dynamics of the dipolar BEC soliton 

[2a - a dd g{K)] + 2£ trap w p , (4) 

[a - a dd h(K)] , (5) 



1 

— + 


1 


N 


W P 


V2^ 




1 

— + 


1 


IN 




V2^ 


w 2 w 2 



with g( K ) = [2- Ik 2 - 4k 4 + 9k 4 cI(k)]/(1 - k 2 ) 2 , h(jc) = [1 + \0k 2 - 2k 4 - 9k 2 cI(k)]/(1 - k 2 ) 2 . The widths of a stationary 
dipolar BEC soliton of energy E = E^ n + E tmp + E{ nt and chemical potential \i = E^ n + E tmp + 2E- mt are obtained by 
solving Eqs. ^ and ^ for w p = w z = 0. 



3. Numerical Results 

We perform numerical simulation of the 3D GP equation ([T]) using the split-step Crank-Nicolson method |[T7l . In 
the presence of the dipolar term the GP equation is integro-differential involving partial derivatives. The dipolar term 
is treated by fast Fourier transformation lfT6l . The error of the reported numerical results is less than 1 %. 

First, we present the results for a single ID dipolar BEC soliton on the square lattice. In this case we take, for 
the weak 2D square OL, Vb = 2. We present the 3D contour of density |0(r)| 2 of the BEC soliton in Fig. [I] for 
(a) g = -3, g dd = 0, (b) g = 0, g dd = 3, (c) g = 1, g dd = 3, and (d) g = 30, g dd = 10 with the 2D square OL 
Vql = _ 2 cos(2x) - 2 cos(2y), where g = 4nNa is the contact-interaction nonlinearity and g dd = 3Na dd the dipolar 
nonlinearity. The corresponding values of number of 52 Cr atoms (a dd = 15ao with ao the Bohr radius (UEl) and the 
scattering lengths are given in the figure caption. The density profiles are distinct in the four cases according to the 
net attraction in the system. The BEC soliton of Fig. [l](b) with no atomic repulsion (g = 0) is the most attractive of 
the four corresponding to a small size, whereas the BEC soliton of Fig Jl] (d) with largest atomic repulsion (g = 30) 
is the least attractive of the four corresponding to a large size. In Fig [lj(a), the BEC soliton is stabilized solely by 
atomic attraction (g dd = 0), and in Fig. [T](b), the stabilization is achieved solely by dipolar interaction (g = 0). The 
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Table 1: Numerical (n) and variational (v) energy and rms sizes, normal-mode frequencies E, (x), (y), (z), Cl p and Q z , respectively, of dipolar BEC 
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Figure 2: (Color online) (a) The phase plot of gdd versus g from variational analysis showing the region of stable dipolar BEC soliton formation on 
square OL with Vo = 2. The *'s denote the numerical points showing the stable-collapse boundary, (b) The numerical (n) and variational (v) rms 
sizes and energy versus g for gdd corresponding to the *'s in (a). 



numerical energy and root-mean-square (rms) sizes of the BEC solitons of Fig. [I] are shown in Table 1 together with 
the variational results. 

One can have a ID nondipolar (gdd = 0) BEC soliton on the 2D square OL for > g > -g C rit, where the numerical 
estimate 8.60 of g cr i t should be contrasted with the variational result of 6.16. As g = 4nNa, these correspond to the 
critical values |Afa|variationai = 0.490 and |Afa|numericai = 0.684. The present numerical critical value is surprisingly 
close to the following critical value when the 2D OL trap in the x-y plane is replaced by the harmonic trap V = p 2 12: 
\Na\ 

numerical — 0.676 1181 . This shows the similar nature of the two nondipolar BEC solitons. 

For a fixed nondipolar nonlinearity g, an ID dipolar BEC soliton (gdd > 0) on the 2D square OL can be stabilized 
for g c Jj tl > gdd > g c Jid 2 > where for gdd < g^ 2 mere is not sufficient net attraction and the BEC soliton expands to 
infinity and for gdd > g^f 1 there is too much net attraction leading to collapse. The domain of stable soliton in this 
case is shown in the gdd versus g phase plot in Fig. [2] (a), where the two lines are the variational boundaries between 
stable BEC soliton and collapse and that between stable BEC soliton and expansion. Of these two lines, the lower 
boundary between stable BEC soliton and expansion can be analytically obtained from the variational equations (|4]) 
and (J5J for cb z - cb p = 0. In this limit the BECs will be infinitely large accommodating an infinite number N of atoms 
and in order that Eqs. ^ and §5§ yield finite numbers for N — > oo one must have [a - addh(K)] = [a - = 0, 

so that H(k) = g(K)/2 with the solution k = 0, while H(k = 0) = 1. Consequently, this boundary is defined by a = add 
or gdd = 3g/(47r) corresponding to the lower line in Fig. [2](a). In this figure the *'s denote the numerically calculated 
boundary between collapse and stability. In Fig. [2](b) we plot the numerical and variational sizes and energies versus 
g for gdd corresponding to the *'s in Fig. [2] (a). 
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Figure 3: (Color online) Numerical 3D contour of the dipolar BEC solitons for (a) g = -3 and gdd = 0, (b) g = and gdd = 3, (c) g = 3 and 
gdd = 3, and (d) g = 15 and gdd = 7 on the triangular OL V^d - _ CO s(2x) - cos(x + V3j) - cos(x - V3y). The density |0(r)| 2 on the contour is 
0.001. 



To perform a linear stability analysis |9] of the BEC solitons on the square OL and obtain the normal-mode 
frequencies, we note that Eqs. ^ and §5§ can be rewritten as ifTUll 

w p = - — , w z = --— , (6) 

ow p ow z 

where U is the linearized effective potential. The squares of the normal-mode frequencies £l z and £l p along axial z 
and transverse directions, respectively, are the eigenvalues of the eigenf unction-eigenvalue problem for the Hessian 
matrix A; 7 = d 2 U/(dwidwj)\ Wz=w ^ Wp=w * where i,j = z,p, where w* and w* are the stationary solutions of Eqs. ^ 
and ([5]) obtained by setting w p = w z = ifTOll . If these frequencies are real, stable oscillation of the widths is 
assured, whereas imaginary or complex frequencies imply exponential increase or decrease of the widths upon small 
perturbation corresponding to unstable states. These frequencies for the states exhibited in Fig. [I] are shown in Table 
1. The real normal-mode frequencies of the BEC solitons of Fig. [T] displayed in Table 1, guarantee their stability. The 
numerical frequencies in Table 1 were calculated from the small oscillations of the widths in real-time propagation 
for a long time. The variational analysis provides a qualitative understanding of many features of the BEC soliton 
including its stability and the normal-mode frequencies. 

Next, we present the results for a ID dipolar BEC soliton on the triangular OL in brief. We take, for the weak 
2D triangular OL, V = 1. We present the 3D contour of density |0(r)| 2 of the stable BEC soliton in Fig. [5] for (a) 
g = -3, gdd = 0, (b) g = 0, gdd = 3, (c) g = 3, gdd = 3, and (d) g = 15, g dd = 7 with the 2D triangular OL 
Vql = _ cos(2x) - cos(x + V3v) - cos(x - V3y). The numerical energy and rms sizes of the solitons of Fig. N 
are shown in Table 2 together with the variational results. The normal-mode frequencies calculated from the linear 
stability analysis |9|, as shown in Table 2, indicate the stability of the solitons. 

The domain of stable dipolar BEC solitons on triangular OL is illustrated in Fig. [4] (a) in a phase plot of the 
nonlinearities g and gdd as obtained from the variational equations ^ and ([5]). Again there is a domain of stable soliton 
between a domain of collapse and of expansion. The collapse takes place for too large a value of dipolar nonlinearity 
gdd and expansion for too small a value of dipolar nonlinearity. A moderate value of the dipolar nonlinearity leads 



Table 2: Numerical (n) and variational (v) energy and rms sizes, normal-mode frequencies E, (x), (y), (z), £l p and Q z , respectively, of dipolar BEC 
solitons of Fig.[3]on triangular OL. 
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to stable solitons. The boundary between stable soliton and expansion is again given by the analytic formula gdd = 
3g/(4n). In Fig. |4](b) we compare the variational and numerical chemical potential and rms sizes (x) and (z) along 
the line gdd — 4 + g/4 covering the whole domain of soliton formation 400 > g > shown in Fig. [4] (a). 

One interesting aspect of studying ID BEC soliton on 2D lattice is to consider an array of many tiny ID BEC 
soliton droplets distributed in 2D OL sites so that a periodic distribution of matter simulating a 2D solid in condensed- 
matter physics with long-range inter-site interaction is obtained. To achieve this, we consider tiny dipolar BEC solitons 
on square OL with g = 1 and gdd = 3 and distribute these on different sites and study the stability of such an array 
by solving the GP equation by real-time propagation. We find that such an array is always unstable due to long-range 
inter-site interaction, if two solitons are placed on neighboring sites along one of the axes - (x,y). However, the 
array is stable if they are placed along the diagonal directions. With this information, we see that a stable periodic 2D 
pattern of tiny solitons is obtained if the occupation of neighboring sites is avoided. At the maximum of a filling factor 
of 1/2, the stable checkerboard configuration is displayed in Fig. [5] (a), where the soliton droplets are put diagonally 
on the black or white spots of a chess board. After real-time evolution of the GP equation, the final array at t = 100 is 
displayed in Fig. [5](d). It is interesting to note that this checkerboard pattern is a stable Mott insulator state obtained 
by solving the Bose-Hubbard model on a 2D lattice with repulsive long-range dipolar interaction ifTTTl . The other Mott 
states obtained there are the stripe and star configurations with 1/3 and 1/4 filling of sites, respectively. We considered 
such configurations with the present tiny dipolar BEC soliton droplets. By reducing the filling factor from 1/2 to 1/3 
or 1/4 we have a lower occupation of sites and hence a lower inter- site interaction. As such arrays get destroyed due 
to the long-range inter-site interaction, a reduced inter-site interaction implies that the stripe and star configurations 




Figure 4: (Color online) (a) The phase plot of gdd versus g from variational analysis showing the region of stable dipolar BEC soliton on triangular 
OL for Vq = I. (b) The numerical (n) and variational (v) rms sizes and chemical potential versus g corresponding to the line gdd = 4 + g/4 in (a). 
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Figure 5: (Color online) Stable array of tiny dipolar BEC solitons, each with g = 1, gdd = 3 on the square OL with Vo = 2 in the (a) checkerboard, 
(b) stripe, and (c) star configurations with filling factors of 1/2, 1/3, and 1/4, respectively at t = 0. The same arrays after real-time propagation at 
t = 100 are shown in (d), (e), and (f), respectively. 



of the BEC soliton droplets are stable. We established the stability of the stripe and star configurations using real-time 
propagation of the GP equation. In Figs. [5](b) and (c) we show the initial stripe and star configurations of the tiny 
BEC droplets at t = 0, and in Figs. [5](e) and (f) the same after real-time propagation at t = 100. The initial and final 
configurations are practically indistinguishable, demonstrating the stability. 



4. Summary and conclusion 

To summarize, using the mean-field 3D GP equation we demonstrated a stable dipolar BEC soliton polarized 
along the axial z direction on a weak square or triangular 2D OL in the orthogonal x-y plane. We considered identical 
OL strengths along different directions and considered a Lagrangian variational analysis of the GP equation with a 
Gaussian ansatz in addition to the numerical solution of the same using the split-step Crank-Nicolson method. The 
stabilization of the BEC was established by the linear stability analysis 0. The widths and energies obtained from 
the numerical solution of the GP equation are in agreement with the corresponding variational results. 

We also considered stable 2D arrays formed by arranging identical tiny dipolar BEC solitons on different sites of 
a 2D square OL. Such an array is unstable if any two adjacent OL sites are occupied. The simplest stable periodic 2D 
array of superfluid dipolar BEC solitons, known as the checkerboard configuration, emerges if the sites are arranged 
diagonally with a filling factor of 1/2. Similar 2D arrays, known as stripe and star configurations for filling factors 
1/3 and 1/4, respectively are also found to be stable. Previously, in the study of ultra-cold dipolar atoms on strict 2D 
lattice using the Monte Carlo simulation of the 2D Bose-Hubbard model, such stable arrays of Mott insulator states 
emerged at filling factors of 1/2, 1/3, and 1/4 [11]. Similar stable configurations obtained in the mean-field GP and 
field-theoretic Bose-Hubbard approaches possibly indicate the general stability property of such structures under the 
long-range repulsive dipolar interaction in the x-y plane. 
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